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We study the configurational properties of single polymers in a theta solvent by Monte Carlo 
simulation of the bond fluctuation model. The intramolecular structure factor at the theta point is 
found to be distinctively different from that of the ideal chain. The structure factor shows a hump 
around q ~ 5/R g and a dip around q ~ 10/ R g in the Kratky plot with R g being the radius of 
gyration. This feature is apparently similar to that in a melt. The theoretical expression by the 
simple perturbation expansion to the first order in terms of the Mayer function can be fitted to the 
obtained structure factor quite well, but the second virial coefficient cannot be set to zero. 
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I. INTRODUCTION 

One of the basic premises in the polymer physics is 
that the bonds connecting neighboring monomer units 
are uncorrelated beyond the persistent length along the 
chain [1]. This property allows us to consider only a 
flexible chain as long as we are interested in large scale 
properties of a polymer chain that is much longer than 
its persistent length. In this sense, it was a little em- 
barrassing to realize that there is actually a long range 
correlation in the bond orientation of a polymer chain, 
and that the bond-bond correlation function decays not 
exponentially but as the power law [2-4] . 

The traditional picture for the bond-bond correlation 
is based on a simple calculation for the polymer with a 
fixed bond angle around freely rotating bonds. In this 
case, one can calculate the bond-bond correlation func- 
tion explicitly to show the exponential decay with a per- 
sistent length [5] . The existence of the long range corre- 
lation, however, has been pointed out [2, 3], and it was 
demonstrated recently that the power law behavior is in- 
duced in the bond-bond correlation through the interac- 
tion between monomers separated by a long distance in 
the curvilinear coordinate along the chain [4, 6-9]. The 
power law in the bond-bond correlation holds not only 
for an excluded volume chain, but also for a chain in a 
melt and in a theta solvent, where a polymer chain is 
supposed to behave as an ideal chain [1, 5]. This has 
been confirmed both by numerical simulations and by 
theoretical analyses. 

This deviation from the ideality of a chain in a melt 
has been seen also in the intramolecular structure factor. 
For the ideal chain, the structure factor decays as q~ 2 in 
the intermediate range, 1/R g < <7 1/a, where R g and a 
are the gyration radius and the bond length, respectively. 
This q~ 2 decay comes from the fractal dimension of the 
ideal chain configuration. The intramolecular structure 
factor of the polymer chain in a melt has been studied 



numerically and theoretically, and it has been found that 
there exists a substantial deviation from the ideal chain 
behavior [8]. 

In this paper, we study the structure factor of a single 
polymer molecule in a theta solvent, that is another situ- 
ation where a polymer chain is supposed to become ideal 
effectively. In a melt, the interaction between monomers 
is screened by the existence of other polymer chains, and 
the excluded volume effect is canceled exactly by the in- 
duced attraction due to the incompressibility of the sys- 
tem [5, 10] while the interaction in a theta solvent is being 
adjusted by some fine tuned external parameter, such as 
temperature, so that the excluded volume effect is com- 
pensated by the attractive part of the interaction. We 
study how this fine tuning of the parameter may affect 
on the virtual ideality of the structure factor. 

This paper is organized as follow. After quickly review- 
ing how the long range correlation comes into the bond- 
bond correlation in Sec. 2, the model and the method of 
our simulations are described in Sec. 3, and the simula- 
tions results are given in Sec. 4. The theoretical analysis 
is outlined in Sec. 5 and the results are discussed in con- 
nection with those for melt in Sec. 6. Detailed expressions 
of the theoretical analysis are given in Appendix. 



II. BOND-BOND CORRELATIONS IN A 
POLYMER CHAIN 

Let us quickly review how the long range correlation 
should arise in the bond-bond correlation along a poly- 
mer chain [4] . We consider a single polymer that consists 
of N monomers. Let r n (n = 1, • • • , N) be the position of 
the n-th monomer in it, and the bond vector is denoted 
as 

a n = r n+ i-r n . (1) 
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We define a subchain as a part of the chain, and introduce 
the subchain vector as 



+ 1 



(2) 



Now, we assume that the bond-bond correlation de- 
pends only on the chemical distance between the bonds: 



P{s) = \(a n - a n+s ) , 



(3) 



where a is the average bond length and (• • • ) denotes the 
ensemble average. Then, the size of the subchain does 
not depend on n, and we have 

R{sf = (Rn(s) 2 ) « a 2 (s + J' dr(s - r)P(r)) (4) 



in the large s approximation. This gives 



(5) 



with lu = 2 — 2v if the subchain size scales as 

R(s) ~ s v (6) 

with the exponent v ^ 1/2. We have lu w 0.824 for the 
excluded volume chain, where v ~ 0.588. 

For the ideal chain with v = 1/2, Eq.(5) gives P(s) = 
as it should, but for the case of apparent ideality of a 
polymer chain in a melt or a theta solvent, we have 

P(s) ~ s- 3/2 (7) 

because of the correction term, 

R(s) 2 « a 2 S (c + c lS - 1 / 2 + ...) . (8) 

III. MODEL AND SIMULATION METHOD 

We perform Monte Carlo simulations of the bond fluc- 
tuation model (BFM) on the three dimensional cubic lat- 
tice [11, 12]. A polymer chain consists of N monomers, 
and the i-th monomer is located at the center of a cubic 
cell T"i, occupying the cell with 8 vertices at lattice sites. 
The bond length between consecutive monomers along 
the chain should be in the range [2, \/l0] with the excep- 
tion of \/8- The i-th. and j-th monomers that are not 
consecutive along the chain, namely j ^ i ± 1, interact 
each other through the "quasi-Lennard-Jones" potential 
energy [13] 
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where = [r, — Tj\ is the distance between the interact- 
ing monomers, and k B and T are the Boltzmann constant 
and the temperature, respectively. The dimensionless pa- 
rameter P is proportional to the inverse temperature and 
characterizes the interaction [14]. 

We perform Monte Carlo simulations, using Metropo- 
lis method along with the slithering snake algorithm to 
accelerate the relaxation towards equilibrium [15]. In our 
system, there is only one polymer chain with N monomer 
units. One Monte Carlo step consists of N trials of ran- 
dom displacement to one of the nearby sites for randomly 
chosen monomers, followed by N slithering snake trials. 

The quantities we study are the radius of gyration R gi 
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the intrachain structure factor S{q) 
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(ii) 
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and the bond-bond correlation P(s) averaged over the 
chain, 



P(s) = 



a 2 (N 
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with the average bond length a, 
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where the angular brackets represent the ensemble aver- 
age. 

Note that the radius of gyration can be expressed us- 
ing the bond-bond correlation by a similar equation as 
Eq.(4): 



1 + 
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P(r) 



(14) 



IV. RESULTS 



A. Theta Point 



First, we have to determined the theta point for our 
model. It is often defined as the point where the second 
virial coefficient vanishes. In our simulations, we define 
the theta point as the point where the radius of gyration 



R g behaves as an ideal chain, 



R g (x VN, 



(15) 



in the large N limit. The interaction parameter (3 C at the 
theta point is determined numerically from the data of 
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FIG. 1: (Color online) The scaling plot for the radius of gy- 
ration R g . The scaled radius of gyration R 2 6/(N"£) 2 are 
plotted against (p - /? C )TV with p c = 0.63, v = <j> = 0.5, and 
£ = 3.17 for various values of TV The inset shows the original 
plot of R 2 g /TV vs p. 
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FIG. 2: (Color online) Bond-bond correlation. P(s)s 3 ^ 2 are 
plotted against s for TV = 1024 for various j3. 



R g for various values of (3 and chain length TV by fitting 
them to the finite size scaling form 



R 2 = In 2 ** f (ip - p D )N+) 



(16) 



with the two exponents, v and <j). The function f(x) 
is the scaling function that satisfies /(0) = 1, and £ is 
the length scale proportional to the average bond length. 
Fig.l shows that the data fit well to the scaling form (16) 
with 

/3 C = 0.63, v = <f> =1/2, ^ = 3.17. (17) 



B. Theta Point determined by Bond-Bond 
Correlation 
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As we have discussed, the bond-bond correlation func- 
tion decays as s~ 3 / 2 at the theta point. Fig. 2 shows the 
bond-bond correlation function P(s) for various temper- 
atures near the theta point for TV = 1024. It shows that 
the data fit to the form 

P(s) s 3 / 2 = B (l- /3//3 c (iV)) 8 + A+-- (18) 

that is consistent with the theoretical results by Shir- 
vanyants, et al. [9]. The theta point may be defined as 
f3 = 0c(N) where the bond-bond correlation decays as 
s~ 3 / 2 , but it turns out that (3 C (N) depends on TV. The 
TV-dependence of (3 C (N) plotted in Fig. 3 shows 

(3 C {N) - p c cx J= (19) 
v TV 

with the theta point in the infinite TV limit f3 c = 0.625. 
This value is close enough to the previous estimate of 
P c = 0.63 by the scaling plot of R g . In the rest of the 



FIG. 3: (Color online) The theta point determined by the 
bond-bond correlation. The values of P where the bond- 
bond correlation shows the ideal behavior, s -3 / 2 , are plotted 
against 1/TV 1//2 . The linear extrapolation to the infinite TV 
shows the theta point, p c — 0.625. 



paper, we use p c = 0.63 for the theta point. Note that 
the TV-dependence of /3 C (N) given by Eq.(19) is consistent 
with the ideal chain behavior of R g of Eq.(15) at = j3 c 
in respect to Eq.(14). 



C. Structure Factor 

The structure factor for the ideal chain, So(q), in the 
large TV limit is given by 

S (q) = Nf D (qR gQ ) (20) 
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FIG. 4: (Color online) Kratky plots for the structure factor 
against Q = qR g . The simulation data at /3 = 0.63 for the 
chain with the length N =1024, 2048, and 4096 are plotted 
along with the Debye scattering function (the dashed line). 
For the value of R g , the values obtained by the simulations 
are used for the simulation data and R g o = Na 2 /6 for the 
analytical expressions. The red curve represent the theoretical 
estimate, Eq.(29), with VNB = -0.55 and A = 0.154. 



with the radius of gyration for the ideal chain 

R% EE \ N a\ 

and the Debye scattering function, 



fo(x) 



1 + x' 



(21) 
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In the intermediate length scale, i.e., l/R g o <C q <C 
1/a, this decays as 



(23) 



The 1 / g 2 -dependence comes from the scaling behavior of 
the ideal chain, therefore, the existence of plateau in the 
plot of S(q) q 2 , i.e., Kratky plot, has been considered to 
be a sign of the ideality in a polymer chain behavior. 

Fig. 4 shows the Kratky plot of our numerical simula- 
tions at the theta point for TV = 1024 ~ 4096; The wave 
number q is scaled by the numerically obtained radius of 
gyration R g as Q = qR g . As the number of monomers N 
increases, the curve in the smaller Q regime tends to con- 
verge to a common trend, but it is clearly different from 
that of the ideal chain (the dashed line). The Kratky plot 
at the theta point shows a hump around Q ~ 5 and a dip 
around Q ~ 10, and does not show the plateau as the one 
would expect for the ideal chain. Their general features 
are apparently similar to those found for a polymer chain 
in a melt [81. 
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FIG. 5: Diagrams that represent contribution for Eq.(28). 



V. THEORY 

Now, we estimate the structure factor theoretically and 
compare it with those obtained by the simulations. Sup- 
pose u(r s ,r) be the interaction potential between the s-th 
and the r-th monomers with r St , r = r s — r r . Then the 
structure factor Eq.(ll) is written as 



S(q) 
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where 
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) represents the statistical average for the ideal 



and (• 
chain. 

We now employ the approximation 



-U/k B T 



using the Mayer function 



f( r r,. 
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Then up to the first order in /, the correction in the 
structure factor SS(q) is given by 

5S(q) EE S(q)-S Q (q) 

2 E E^^/k^o 
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iq-r„,„ 



> Q (/(r r ,.)> • (28) 



There are four types of contribution in Eq.(28) depend- 
ing upon the relative positions of n, m, r, and s (Fig. 5). 
Adopting the bead-spring model for the ideal chain aver- 
age (...) with the average bond length a, and employing 
the further approximation valid for ga C 1, we obtain 
the expression 



N 



SS(q) w F(qR g0 )VNB + G(qR g0 )A (29) 
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FIG. 7: (Color online) The large Q behavior of the structure 
factor after subtracting C/Q 2 with C = 1.336. 



clearly that the second leading order is 1/Q 3 , which is 
consistent with the asymptotic expression we obtained 
in Eq.(A20). 



FIG. 6: (Color online) Theoretical results for structure fac- 
tors. The numerical estimates by Eq.(29) are plotted with 
A = 0.154 and VNB = -0.55, and 0. 

up to the second leading order in N. The functional 
forms for F and G are given in Appendix. The dimen- 
sionless parameters, 

B = ™ J drf(r), (30) 

A ee 1 J drr 2 f(r), (31) 

characterize the interaction. The parameter B is twice 
of the second virial coefficient for the unlinked monomer 
gas [16], and is supposed to be close to zero for the theta 
solvent [10]. Note that B comes into S(q) as \/~NB. 

Fig. 6(a) shows the Kratky plots for the theoretical 
structure factor S(Q)Q 2 with Q = qR g o by Eq.(29) with 
A = 0.154 and ^/~NB = —0.55 along with the curve for 
A = 0.154 and B = 0, the Debye function, and the simu- 
lation data for N = 4096. One can see that the curve for 
B = is almost proportional to the Debye function and 
cannot be fitted to the simulation data. Fig. 6(b) shows 
SS(Q)Q 2 /N with the simulation data for various values 
of N. The data converge to the theoretical curve quite 
well as N increases. 

In Fig. 7, the large Q behavior of the structure factor 
are plotted in the logarithmic scale, after subtracting the 
leading order term of C/Q 2 with C = 1.336. It shows 



VI. DISCUSSIONS 

Our findings are summarized as follows: (i) The theta 
point for the finite chain (3 C (N), where the bond-bond 
correlation decays as s~ 3 / 2 , depends on the chain length 
A, and in the infinite chain length limit it converges to 
the theta point j3 c determined by the scaling behavior of 
the radius of gyration, (ii) The structure factors at the 
theta point (3 C obtained by the numerical simulations are 
distinctively different from that of the ideal chain, i.e., 
(ii-a) in the intermediate range of q > 1 / R g , Kratky plot 
of the structure factor shows a hump and a dip, and (ii- 
b) in the larger q range, 1/R g -C q -C 1/a, the structure 
factor decays as 1/q 2 with the next order term of 1/q 3 - 
(iii) In the large N limit, numerically obtained structure 
factors fit well to the simple perturbation expression up 
to the first order of the Mayer function with the fitting 
parameters VNB = -0.55 and A = 0.154. 

In comparison with the case of a melt, the obtained 
structure factors at the theta point apparently resemble 
the ones of a melt in the existence of a hump and a dip 
[8]. In the large q region, however, the structure fac- 
tors for the two cases differ; In the case of melt, it has 
been shown by numerical simulations that structure fac- 
tors decays as 1/g 3 for the large q region [8], in contrast 
with the present case, where we found the 1/q 2 decay 
with the positive 1/q 3 correction. The 1/q 3 decay in the 
melt was interpreted as a result of renormalization from 
the 1/q term obtained by the one-loop approximation [8] 
while the 1/q 2 term and the 1 /q 3 term at the theta point 
directly correspond to the A-teim and the V^/VS-term, 
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respectively, in our theoretical expression. In the real 
space, the 1/q 2 decay of the structure factor in large q 
means the 1/r density correlation of the ideal chain in 
short r. The existence of the positive 1/q 3 correction 
and the hump around q ~ 5/R g implies that the density 
correlation does not decay as fast as that of the ideal 
chain around r ~ R g . 

The simulated structure factors here at the theta point 
resemble those of a melt, but it is intriguing that we 
need to set y/NB ^= for our theoretical expression of 
the structure factor to fit to the numerical results at the 
theta point [3 C . The theta point is often considered as the 
point where the second virial coefficient B/2 vanishes. 
Actually, if we set B = 0, the structure factor shows 
Kratky plateau and looks pretty much like the ideal one 
(the dotted line in Fig. 6(a)), but this cannot be fitted 
to the simulation results by adjusting A alone even with 
an arbitrary factor. Theoretically, if we consider higher 
order corrections, the theta point should correspond to 
the vanishing point of the interaction parameter z 

z = VNb (32) 

where b is the excluded volume parameter with the ad- 
ditive renormalization by the higher order cluster contri- 
butions, and may be given in the form 

b = B + CN- 1 ' 2 + ... (33) 

with a constant C . However, such renormalization effect 
could be partially taken into account by replacing B with 
b, and one would expect the structure factor at the theta 



point should be given by the one with ^/NB — in our 
expression. 

In the bond-bond correlation P(s), we observed the 
analogous ./V-dependence, which can be interpreted in 
the same way. Shirvanyants et al. [9] have obtained the 
theoretical expression corresponding to Eq.(18) in a sim- 
ilar approximation to the one we employed for S(q), and 
showed that the coefficients for the s -1 / 2 and the s~ 3//2 
terms in Eq.(18), i.e. B (l - (3/ (3 C (N)) and A, are pro- 
portional to the parameters B and A, respectively. The 
observed TV-dependence in (3 C {N) as Eq.(19) in our sim- 
ulations should be a result of higher order effects, that 
could be obtained by replacing the coefficient of s~ 1//2 - 
term with the renormalized one as Eq.(33). 

Regarding the origin for the deviation of the structure 
factor from the Debye function, an obvious possibility 
could be that inaccurate estimate for the theta point (3 C , 
i.e. our estimate f3 c — 0.63 is not close enough to the 
theta point for the structure factor to be Debye-like even 
though the two independent estimates from the radius 
of gyration R g and the bond correlation coincide with 
each other within our numerical precision. This might 
happen due to the slow convergence caused by the loga- 
rithmic correction at the theta point. Another possibility, 
perhaps more interesting one, would be that the tricriti- 
cal fluctuations at the theta point produce a non-trivial 
contribution to the structure factor. 

Acknowledgments 

The authors thank Dr. Takahiro Sakaue for critical 
reading of the manuscript. 



APPENDIX A: FIRST ODER CALCULATION OF S(q) IN THE MAYER FUNCTION 



In the appendix, we describe the calculation of Eq.(28) to obtain the explicit expressions for F and G in Eq.(29). 
The contribution from the diagram (i) in Fig. 5 is given by 



Si(Q) = ± E [(^• r -/K s )) -(e^-) (/(r r;S )) 
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where 



H(q;h,l 2 ) = J dn/(n) J dr 2 e^ r2 (Go(r 1 -r 2 ;/ 1 -? 2 )-Go(r 1 ;? 1 ))Go(r 2 ;/ 2 ) 
with the free propagator 
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H(q; li,h) can be estimated as 
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In the last expression, we have expanded the cosine up to q 2 7 which is valid for qa -C 1. 
The contributions from (ii), (hi), and (iv) can be given in terms of H as 

Su (q) = A £ e -(l/6), a a a (n-r) ff(g . r _ ajr _ s)e -(l/6), a a a (.-m) 



l<m<s<r<n<<JV 



r ^H(q; r — s,r — m) 



Siii{q) = e ~ 

l<s<m<r<n<<N 

= S iv (q). 

For large N, we can replace the summation by integral and obtain, to the leading orders in N, 

Si(q) = N 3 / 2 Fi(qR g0 )B + N 1 / 2 G l (qR g v)A, 
S u {q) = N 3 ^ 2 F il (qR g ())B + NG il (qR g0 )A, 
(q) = N 3 / 2 F ui (qR g0 )B + N^ 2 G iu (qR g0 )A, 
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where C(s) is Riemann's zeta function and C(3/2) = 2.61237- • •. 

With these functions, the correction of the structure function is given by 
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5S(g) « F(qR g0 ) VNB + G(qR g0 )A 
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with 



F(Q) = F i (Q) + F ii (Q) + 2F iii (Q), 
G(Q) ee Gn(Q) 

in the leading orders in N. For large Q 7 i.e., l/R g o <C q <C 1/a, we have 
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with C w 1.01171, and for small Q, i.e., q < 1/i? 
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From Eq.(A21), the correction for the radius of gyration R g is obtained as 
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